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ABSTRACT 

We analyze the light curves of 413 radio sources at submillimeter wavelengths 
using data from the Submillimeter Array calibrator database. The database in¬ 
cludes more than 20,000 observations at 1.3 and 0.8 mm that span 13 years. We 
model the light curves as a damped random walk and determine a characteristic 
time scale r at which the variability amplitude saturates. For the vast majority 
of sources, primarily blazars and BL Lac objects, we find only lower limits on 
r. For two nearby low luminosity active galactic nuclei, M81 and M87, however, 
we measure r = days and r = 45^24 days, respectively (2 a errors). In¬ 

cluding the previously measured r = 0.33 ± 0.16 days for Sgr A*, we show an 
approximately linear correlation between r and black hole mass for these nearby 
LLAGN. Other LLAGN with spectra that peak in the subnnn are expected to 
follow this correlation. These characteristic time scales are comparable to the 
minimum time scale for emission processes close to an event horizon, and suggest 
that the underlying physics may be independent of black hole mass, accretion 
rate, and jet luminosity. 
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Introduction 


In recent years there has been much work focused on understanding black hole accretion 
and its similarities across the mass scale, from stellar mass black holes in X-ray binaries 
(BHBs) to active galactic nuclei (AGN). Long term variability studies have found evidence 
for a mass-dependence in timing features that holds from BHB to AGN (McHardy et al. 2006 


Kelly et al.|2009 ; MacLeod et al.|2010 ). A linear correlation can be understood if the emission 
originates in a region that is the same size for all systems in units of Schwarzschild radii, 
where Rs = 2 GM/c 2 . However, the timescale is also inversely proportional to the Eddington 
accretion fraction, suggesting that the pertinent size of the emission zone is also regulated 
by the total system power. Similarly, studies of broadband spectra have found a strong 
correlation between radio and X-ray luminosity and black hole mass (the “Fundamental 
Plane of Black Hole Accretion”, or FP, see e.g., |Merloni et al.|[2003| |Falcke et al.||2004 


Plotkin et al.||2012 ), that have confirmed earlier theoretical frameworks (Blandford & Konigl 


1979 Falcke & Biermann 1995) that synchrotron spectral features scale predictably with 


black hole mass and accretion power. Combining these concepts, observations at a fixed 
frequency of black holes of the same mass with different accretion power should not yield 
similar timescales, because the frequency “selects” out different sized emission regions in 
both cases. 

The sub-millimeter (submni) band seems to be selecting out regions of event-horizon 


scale in two sources of drastically different mass and accretion rate: Sgr A* and M87 (Doelc- 


man et al. 2008, 

2012afb 

). Both sources belong to the class of nearby, low-luminosity AGN 

(LLAGN; Ho 

1999 

) that fall on the FP. 


Variability provides an independent test of the size scales. Dexter et al. (2014) used 


submm light curves to demonstrate that Sgr A* follows a damped-random walk (DRW) 
variability pattern, with a characteristic time scale r = h at 230 GHz at 95 per cent 
confidence, with consistent results at higher frequencies. This time scale is an order of 
magnitude larger than the period of the last stable orbit for a non-rotating black hole, which 
is most easily understood as resulting from accretion processes on a scale of a few to 10 R s . 

We present a variability study of AGN in the submm, including the LLAGN M81 and 
M87 as well as the 411 other radio sources included in the SMA Calibrator database, which 
span more than a decade in duration (§21). This sample was previously analyzed by Strom 


et al. (2010), but with a much smaller number of objects and total span. The full sample of 


light curves are well described by a noise process, which we model with a damped random 
walk in order to quantitatively measure characteristic variability time scales (§[3]). We show 
that the LLAGN follow a correlation between black hole mass and time scale, while the 
higher luminosity AGN have much longer characteristic time scales with no clear black hole 
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mass dependence. 


2. Data Sample 


The SMA calibrator database (Gurwell et al. 2007) includes observations of 412 sources 
from June 2002 to January 2015. A total of 23254 flux densities are recorded with 19111 of 
those flux densities obtained in the 1.3mm band. 410 light curves have > 10 flux densities, 
141 light curves have > 30 flux densities, and 48 light curves have > 100 flux densities. 
Observations were obtained at frequencies from 200 to 406 GHz with 70% obtained between 
220 and 235 GHz. 165 observations (147 at 1.3mm) of M87 are included in the calibrator 
database spanning from January 2003 to January 2015. Data reduction of sources followed 
standard SMA techniques that set flux densities on absolute flux density scales determined 
by solar system objects with typical accuracy of 5 to 10%. 

Sources are selected for the SMA calibrator database on the basis of their suitability 
as phase reference sources for mm/submm interferometric observations, primarily reflected 
in their compact size, bright radio flux density, and declination above -50°. Flux densities 
in the data base range from 26 mJy to 52 Jy with a median value of 1.2 Jy. 336 sources 


are matched to flat spectrum radio sources in the CRATES catalog (Healey et al. 2007). 


Non-matches to CRATES include steep spectrum calibrators such as 3C 286 and sources in 
the Galactic plane such as J1700-261 that may have been missed by the single dish surveys 
that form the basis of CRATES. The CRATES sample includes primarily blazars and BL 
Lac sources; 289 of the sources are found in the CGRABS catalog of y-ray blazar candidates 


(Healey et al. 2008). 


We also included SMA monitoring observations of M81 in our analysis, which has been 
the target of a significant campaign (McHardy et al. 2015). This campaign included 86 ob¬ 
servations between September 2009 and March 2012, of which 42 observations were obtained 
at frequencies higher than 300 GHz. We also include PdBl observations of M81 obtained in 
2005 at 1.3 mm with a time resolution of ~ 1 hour (Schodel et al. 2007). 


3. Time Scale Analysis 

The light curve for M87 is shown in Figure [l] The light curve for M81 is shown in 
McHardy et al. (2015) and shares similar characteristics of short-term variability and long¬ 
term stability. These properties are visible in the structure functions (Figj2|, defined for a 
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Fig. 1.— Light curve for M87 at 1.3 mm from the SMA. The light curves show the charac¬ 
teristic of short-term variability coupled with long-term stability that is the signature of the 
damped random walk. 
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light curve s(t) with N data points as 


S 2 (At) = -E( S (t)- S (f + At)) 2 . 


( 1 ) 


Both structure functions show increases in activity on short time scales followed by a 


plateau, similar to those of Sgr A* (Dexter et al. 2014) but markedly different than those of 


typical calibrators in the database. For instance, the structure functions for 3C 84 and 279 
rise smoothly without any apparent saturation time scale for At < lOyr. Structure functions 
are a useful guide for qualitative comparison; however, they functions are unreliable for 
quantitative analysis ( Emmanoulopoulos et al.||2010 ). They can show artifacts from aliasing 
(e.g. peaks at ~ 180 days associated with the anual cycle), as well as spurious saturation time 
scales even with regularly sampled data. Fitting of structure functions for saturation time 
scales is also challenging since their values and errors at different time scales are correlated. 


3.1. Damped random walk model 


In order to quantitatively measure characteristic time scales, we model the entire set 
of subnim light curves as the result of a stochastic damped random walk (DRW) process. 
The DRW is a simple 3-parameter model (mean, standard deviation of long-timescales, and 
transition time) that is well-suited to parametrize the most important properties of noise 


light curves (Kelly et al. 2009). The DRW consists of red-noise on time scales less than r 


and white-noise on longer time scales. The resultant structure function for the DRW is 

S\t) = Si (1 - , 


( 2 ) 


where S ^ is the power in the light curve on long time scales At r. The characteristic 
time scale r determines the de-correlation time on which the variability transitions from red 
noise on short timescales to white noise on long time scales. This model has been shown to 


successfully describe the optical light curves of quasars (Kelly et al. 2009 MacLeod et al. 


2010), as well as the submm light curve of Sgr A* (Dexter et al. 2014). 


We convert the likelihood of the observed light curves arising from a given set of DRW 
parameters into the posterior probability of the set of model parameters using a Bayesian 
approach. Kelly et al. (2009) and Dexter et al.| (2014) used a Metropolis-Hastings algorithm 
to sample the likelihood over the model parameter space. We use this approach to measure 
the model parameters in cases where the upper limit to r is smaller than the light curve 
duration T. However, this is not the case for the vast majority of the SMA calibrator light 
curves. When no upper limit can be found, the Metropolis-Hastings algorithm fails, as it 
preferentially samples the long tail of the probability distribution rather than its peak. 
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To avoid this issue, we sample the probability distribution over a regular 64 3 grid in 
the parameters with uniform priors for the mean (//), standard deviation (log cr), and logr. 
Since for the vast majority of sources we are interested in lower limits on r, we choose the 
parameter ranges to ensure that the lower end and peak of the probability distribution are 
well captured, at the expense of the long tails to large r (much longer than the duration of 
the light curve). From the probability distributions in r, we consider good limits to be cases 
where r > At and r <T at 99.7% confidence, where At is the minimum separation between 
two measured points in the light curve and T is the total light curve duration. In all cases 
we report 2a limits because of the non-Gaussian nature of the distributions that can include 
long tails. 


Although less efficient than Monte Carlo, this direct grid method recovers consistent 
parameter estimates for the M81 and M87 light curves. If the true light curve is not described 
by the DRW, as seen in Kepler quasar light curves (Kasliwal et al. 2015), the resulting r 
estimate will be biased. For light curves that can successfully be fit by the model (xt ~ 1), 
such as the SMA data used here, the bias leads to < lcr changes in our estimates of r. 


3.2. Results 


We use the DRW model to infer r values from the full set of 664 230 and 345 GHz light 
curves from the 413 sources. A histogram of the resulting 95% confidence lower and upper 
limits on r is shown in Figure [3] In total, we find 276 light curves with lower limits (40%) 
and 16 with upper limits (2%). The histograms are further broken down by sources with 
> 30 data points and those with available black hole mass estimates in the literature (see 
below). For the measured lower limits, the sub-samples are consistent with the full sample. 
Many of the measured upper limits are false positives in light curves with few data points 
caused by not sampling a large enough range in r. After those are removed, 5 detections of 
upper limits remain, all of which have xl close to unity, similar to Sgr A*. In Table [lj we 
summarize the detections. The light curve residuals after model subtraction are normally 
distributed for M81 and M87, as expected for an accurate model, and similar to what was 
seen for Sgr A* (Dexter et al. 2014). 


The results of this time scale analysis are in stark contrast for the available LLAGN light 
curves (Sgr A*, M81, M87) and those of ordinary AGN / blazars (the vast majority of the 
light curves). For the LLAGN, 3/3 sources have well measured values of r. For the rest of 
the sources, 1% (4/664) have reliable measurements of r. It is therefore potentially possible 
to select LLAGN based on submm properties alone, with a low false positive rate and a high 
detection efficiency for well sampled light curves. There is no systematic difference in the 
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number or cadence of measured flux densities that can explain this result: it is related to an 
intrinsic difference in the submm variability properties of the different source classes. 


3.3. Black Hole Mass Estimates 


We compiled black hole masses from the literature. The three nearby LLAGN have the 
most accurately determined masses, often from multiple dynamical methods. Sgr A* has a 


mass of 4.4 x 1O 6 M 0 ( 

Genzel et al. 

2010’). 

M87 has a mass of 3.5+q; 9 x 10 9 M 0 from gas 

dynamic measurements ( 

Walsh et a 

2013 

) and a mass of 6.6 ± 0.4 x 1O 9 M 0 from stellar 


dynamic measurements (Gebhardt et al.||201 1). Following Kormendy & Ho (2013), we adopt 


a value of 6.2 x 10 9 M 0 . Similarly for M81, stellar and gas dynamical measurements exist 


(Bower et al. 2000 Devereux et ah 2003) and we follow Kormendy & Ho (2013) to adopt a 


value of 6.5 x 1O 7 M 0 . 

For the remainder of the sample, we rely on published black hole mass estimates from 


various techniques. These include broad line region spectroscopy (Woo & Urry 2002 Kelly 


& Bechtold 2007 

Shen et al. 2008, 

2011 

Paris et al. 2014 

), black hole mass-bulge luminosity 


correlation (Woo & Urry 2002), and the FP between radio-X-ray luminosity (Merloni et al. 


2003 Plotkin et al. 2012). For the case of Paris et al. (2014), which does not provide mass 


estimates, we use their measurements of the Mg II transition and the methodology of Shen 


et al. (2011) to estimate masses. FP estimates of the black hole mass are derived using 
ROSAT 0.2 - 2.0 keV X-ray fluxes and the 1.3mm flux density that is debeamed assuming a 
Doppler boosting factor of 7. We identify 194 black hole mass measurements for 148 sources. 
For sources with multiple mass estimates, we rely on the most recent published estimate. 
The accuracy of the mass estimates range from 0.16 dex for Mg II transitions to > 0.4 dex 


for continuum luminosity (Shen et al. 2011) 


4. Discussion 

The t-Mbh relationship for the sub-sample of SMA calibrator sources with black hole 
mass estimates is shown in Figure |4| The difference in submm variability properties between 
regular and low-luminosity AGN is immediately apparent: the vast majority of AGN only 
yield lower limits on r, while the LLAGN have measured time scales that lie well below 
many of these at comparable black hole mass. 

On the other hand, the measured r values for the 3 LLAGN with good measurements of 
r do show a significant trend of increasing r with increasing black hole mass. Formally, the 


























































Fig. 2.— Structure functions for M87, M81, 3C 279, and 3C 84 from 1.3 mm SMA data. 
The structure function suggests that M87 and M81 have different variability characteristics 
from the two high power radio sources. Quantitative analysis with the DRW model confirms 
the differences. 
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log[x (d)] 



Fig. 3.— Histogram of r lower (left) and upper (right) limits for entire sample (black), 
those with n > 30 (blue), and those with n > 30 and black hole mass measurements (red). 
Constraints with n < 30 are used to avoid false positives. Those with black hole masses are 
representative of overall distribution. 


Table 1. Characteristic Time Scales for Submm Variability 


Source 

r 

"How 

Thigh Xu 

Mbh 

Source 


(days) 

(days) 

(days) 

(Af®) 



0433+053 

106 

63 

299 

0.88 

2.6 x 10 7 

Woo & Urry 

(200 

2) 

0721+713 

75 

56 

155 

0.83 

5.2 x 10 6 

Fundament a 

da lie 

1104+382 

40 

16 

154 

0.81 

1.9 x 10 8 

Woo & Urry 

(200 

2) 

1751+096 

106 

50 

331 

1.04 

2.2 x 10 7 

Fundament a 

da lie 

M87 

45 

21 

106 

0.87 

6.2 x 10 9 

Kormendy & Ho 

2013) 

M81 

1.6 

0.7 

4.6 

1.14 

6.5 x 10 7 

Cappellari et al. ( 

2009) 

Sgr A* 1 

0.33 

0.17 

0.5 

1.05 

4.4 x 10 6 

Genzel et al. 

(201 

0) 


Dexter et al. (2014) 
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Fig. 4.— Black hole mass versus DRW time scale for the SMA calibrator source sample. 
LLAGN Sgr A*, M81, and M87 are shown as filled in circles. Other AGNs from the SMA 
calibrator database with black hole masses are shown as open circles for detections and up 
arrows for lower limits on r. The time scales associated with AGN/blazars are significantly 
longer at the same black hole mass than those of the LLAGN. The inferred time scales for 
LLAGN variability increase with black hole mass, and are consistent with the predicted 
linear relationship from General Relativity. The solid line gives the period for the innermost 
stable circular orbit (ISCO); the dashed line gives the infall time for a disk with radius 5 Rs, 
a viscosity a — 1, and a height H/R = 1. 
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best power-law fit is r ~ 0.3 day (Mbh/4 x 10 6 M o )°’ 7±(U (with a coefficient of determination 
R 2 = 0.996). We consider the observed trend consistent with a linear relationship, however, 
given the small number of sources and the large uncertainties in r and Mbh- If the M87 
mass were taken from gas dynamical measurements, for instance, then a linear law would 
be consistent with the data. The reliability of the correlation rests on only three LLAGN. 
This is a minimally small number of objects for a correlation, but these are the only LLAGN 
for which r can currently be determined. We are performing observational campaigns to 
increase the number of LLAGN with measured r. 

A prediction of general relativity is that time scales around black holes should scale 
as r ~ A 3//2 Mbh, where R is the radius in units of the event horizon size. This time scale 
reaches a minimum for a given Mbh when R is comparable to the size of the event horizon, 
and in that case one might expect to find a linear relationship between r and Mbh- This is 
clearly not borne out by the full sample: there is a large scatter in measured r values for 
regular AGN, and they do not fall on a linear track with the LLAGN. The discovery of such 
a relationship for the LLAGN suggests that we are measuring such minimum time scales 
from emission close to an event horizon. 


This interpretation is consistent with our understanding of these sources inferred from 
their spectral energy distributions. For all three LLAGN, the spectra peak in the millimeter 
or sub millimeter regime ( Bower et al.|2015 Markoff et al. 2008 Doi et al.|2013 ). Synchrotron 
self-absorption leads to a large photosphere which shrinks until the spectral peak, where we 
can see down to the event horizon. At longer wavelengths, we expect longer timescales 
because of the larger scales of the photosphere. This larger photosphere is seen as a growing 
intrinsic image size with wavelength in VLB1 observations of Sgr A* (Bower et al. 2006) and 
M87 (e.g., Hada et al. 2011). The same effect shows up in the variability time scales for Sgr 
A* and M81 and likely for M87. For Sgr A*, Macquart & Bower (2006) found a characteristic 
time scale at wavelengths from 0.7 to 20 cm that ranged from 6 days to hundreds of days, 
increasing with wavelength. Additionally, we performed a DRW analysis of the 2cm M81 
light curve (Ho et al. ||1999 ) and found a lower limit to r of tens of days, much larger than the 
measured value at 1.3mm. In addition, for the case of Sgr A*, the NIR variability time scale 
matches the submm value, consistent with the hypothesis of emission at both wavelengths 
emerging from the same region near the event horizon (Witzel et al. 2012). Thus, one should 
expect a similar relationship for other LLAGN with spectra peaking in the submm. The 
blazars and high power radio sources presumably do not follow a linear correlation because 
these sources are still optically thick at this wavelength. X-ray binaries in the low/hard 
state may also show the same correlation in the optically thin regime, although this may be 
difficult to observe due to the ~ 1 sec time scales implied by the relationship and the rapidly 
evolving dynamics of these syterns (Russell et al. 2014). 
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A similar relationship was found in the X-ray (McHardy et al. 2006) for black holes 


ranging from stellar to supermassive. This relationship, however, relied on a scaling factor 
based on the accretion rate, which we do not require here. This scaling factor could be 
interpreted as compensating between different values of R in the time scale corresponding to 
the X-ray emission region. In the LLAGN, we instead appear to have reached the minimum 
variability time scale corresponding to emission from event horizon scales around black holes. 

The scaling relationship between black hole mass and variability time scale for LLAGN 


is an important insight for Event Horizon Telescope observations of these sources (Fish et al. 


2013). These imaging observations will have resolutions as good as a few Schwarzschild 
radii and have the potential to probe fundamental gravitational physics (e.g., Lu et al. 


2014 Broderick et al. 2014; Pen & Broderick 2014; Ricarte & Dexter 2015). The r — M BH 


correlation shows that the compact sizes measured for these sources are actually colocated 
with the black hole. In addition, many of the fundamental parameters such as spin are 
degenerate with source physics parameters including the accretion rate, the jet inclination 
angle, the electron temperature distribution, the magnetic field strength and orientation, and 
the optical depth. Given the vastly different scales, environments, and physical properties 
of these three LLAGN systems, the existence of this correlation demonstrates a surprising 
coherence among these sources. These, and potentially other LLAGN with spectra peaking 
in the mm/submm can be treated as a unified class. Broadband spectra and light curve 
monitoring of other nearby LLAGN in the submm will add to the number of sources that 
fall into this class. Sources that do not follow this correlation, such as the majority of the 
SMA calibrators, may follow the relation at a higher frequency where the emission becomes 
optically thin. 


The Submillimeter Array is a joint project between the Smithsonian Astrophysical Ob¬ 
servatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded 
by the Smithsonian Institution and the Academia Sinica. This research has made use of the 
NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Lab¬ 
oratory, California Institute of Technology, under contract with the National Aeronautics 
and Space Administration. 
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